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A dynamical method for inelastic transport simulations in nanostructures is compared with a 
steady-state method based on non-equilibrium Green's functions. A simplified form of the dynami- 
cal method produces, in the steady state in the weak-coupling limit, effective self-energies analogous 
to those in the Born Approximation due to electron-phonon coupling. The two methods are then 
compared numerically on a resonant system consisting of a linear trimer weakly embedded between 
metal electrodes. This system exhibits enhanced heating at high biases and long phonon equilibra- 
tion times. Despite the differences in their formulation, the static and dynamical methods capture 
local current-induced heating and inelastic corrections to the current with good agreement over a 
wide range of conditions, except in the limit of very high vibrational excitations, where differences 
begin to emerge. 



I. INTRODUCTION 

The effects of inelastic interactions between current- 
carrying electrons and the vibrational motion of atomic 
nuclei is one of the principal phenomena of interest in 
the field of molecular electronics. These effects have been 
extensively studied experimentally in recent years 
3] . The operation and properties of atomic-scale devices 
are strongly dependent on electron-ion interactions. The 
inelastic scattering of electrons by nuclei, and subsequent 
dynamical motion of the atoms, influences the transport 
properties of the device, while local Joule heating within 
the junction limits the stability of the device. 

The simplest approach to such phenomena is lowest- 
order electron-phonon scattering theory, i.e. the Fermi 
Golden Rule (FGR). This includes the first order correc- 
tions to the electronic system from the electron-phonon 
interaction, which is treated as a perturbation. Phe- 
nomena such as the injection of power in the vibrational 
modes of atomic wires [H and corrections to the current- 
voltage spectrum which arise from the presence of in- 
elastic electron-phonon scattering can be captured at a 
qualitative level within this framework. First-order per- 
turbation theory however cannot be expected to handle 
the limit of strong electron-phonon coupling, or the ef- 
fects of multiple scattering. 

An established method of generalising the FGR to in- 
clude higher-order processes is non-equilibrium Green's 
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function theory (NEGF) g, 0]. One manner in which 
this is done is to consider only lowest-order Feynman dia- 
grams in the expression for the self-energy and to expand 
the Dyson equation in a Born series in the free Green's 
functions. If the electronic Green's function used in the 
Dyson equation and in the calculation of the self-energy 
are the same, one obtains the Self-Consistent Born Ap- 
proximation (SCBA). SCBA has been applied to inelas- 
tic transport both in model systems [^, and, to- 
gether with first-principles electronic-structure calcula- 
tions, in realistic atomic chains and molecular-wire sys- 
tems [ni[il[ll. The SCBA technique is outlined further 
in the following section. The Green's function method 
can be applied also in the time domain, in order to take 
account of transient effects and the response of the sys- 
tem to dynamical driving fields [3, ■ 

Recently, an alternative method for inelastic trans- 
port has been proposed [3, [13, [13 that differs from 
NEGF in philosophy and formulation. The key aim of 
this method is to extend molecular dynamics by rein- 
stating electron-nuclear correlations and the quantum 
nature of nuclei in order to produce a computationally 
tractable form of quantum correlated electron-ion dy- 
namics (CEID) that retains inelastic electron-phonon in- 
teractions, energy transfer and dissipation between the 
two subsystems. Thus far, the method has been ap- 
plied to inelastic I — V spectroscopy in atomic wires 
[l7| and, when combined with electronic open bound- 
aries, was used to calculate local heating in atomic wires, 
and its signature on the current, in real time [l^. An 
outline of the method is given in Section IIIII 

In this paper, we report the first direct comparison of 
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the dynamical CEID method with the SCBA. For this 
comparison, we have chosen a particular model system 
that exhibits interesting behaviour. We study a linear 
trimer weakly coupled to metal electrodes with the cen- 
tral atom in the trimer allowed to move. In the absence 
of electron-electron screening and coupling of vibrations 
to the surrounding lattice, this resonant system is found 
to undergo local Joule heating that is significantly larger 
than that obtained in a ballistic wire. The time taken 
for local phonons to equilibrate with the current-carrying 
electrons is also enhanced, and is strongly dependent on 
voltage. 

The outline of the present work is thus as follows. In 
the next section, the SCBA formalism is outlined. Fol- 
lowing that, the CEID methodology is outlined in Sec- 
tion [TITl and it is shown, to lowest order in the electron- 
phonon coupling, that the steady-state solution to the 
one-particle electronic density matrix involves effective 
sclf-encrgies which are analogous to those in the Born 
Approximation. We also examine the infinite-mass limit 
of the CEID equations, and demonstrate that they re- 
duce to the exact solution of a specific clastic scattering 
problem. The combination of CEID with electronic open 
boundaries is briefly summarized. 

In Section [TVl the static and time-dependent methods 
are applied to our model system. The inelastic I ~ V 
spectrum is analysed using both methods for two limit- 
ing regimes; one where (i) the moving ion is assumed to 
remain always in its ground state (the externally damped 
limit with perfect heat dissipation to the electrodes), and 
the other where (ii) no lattice heat conduction is allowed 
(the externally undamped limit with maximal heating). 
The inelastic current as a function of the thermal excita- 
tion of the quantum ion is studied for a variety of ionic 
masses, and the methods agree up to ionic vibrational 
energies ~ 1 eV. Differences that emerge under more ex- 
treme conditions, and other directions for future work, 
are summarized at the end. 



II. THE SELF-CONSISTENT BORN 
APPROXIMATION (SCBA) 

In this section, the formalism of the SCBA is briefly 
outlined. The detailed description of the method is 
outlined elsewhere d, H, One assumes a coupled 

electron-phonon system within the harmonic approxima- 
tion, whose Hamiltonian, in second quantization, is writ- 
ten as 
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Here Hq is the electronic Hamiltonian described via the 
one-electron basis evaluated at the classical equi- 

librium nuclear positions {Ro}, and {cj^^^} is a corre- 
sponding set of one-electron annihilation (creation) oper- 
ators, -ffph is the phonon Hamiltonian for a set of uncou- 
pled harmonic oscillators with {a^^"*} the set of annihila- 
tion (creation) operators within the occupation number 
representation, and fl\ is the vibrational frequency of 
mode A. H^-ph describes the interaction between the 
electron and phonon subsystems, where the matrix Af ^ 
is the electron-phonon coupling matrix for phonon mode 
A. We also impose the non-crossing approximation, as- 
suming that the interaction of the electron gas with the 
electron reservoirs is independent of its interaction with 
the vibrational modes of the system. 

We further assume that the electron Green's functions 
Gq'-^ for the phonon-free electronic system can be evalu- 
ated. In the case of a nanoscale system coupled to exter- 
nal electronic reservoirs, these will explictly include the 
contribution due to the device-electrode coupling. The 
bare phonon Green's functions Dq'-^ in the frequency do- 
main are those of a free harmonic oscillator of frequency 

^^^""^ = c.-l^o + i^/'c. + fio + W 

Dfiiu) = -2Tri[{Nph + l)6{Lu±no) + Np^Siu;T^om 

where 77 ^ 0"*" and A^ph is the phonon occupation num- 
ber which in equilibrium is given by the Bose-Einstein 
distribution. 

In the weak coupling limit, it is appropriate to con- 
sider only the lowest-order phonon contributions to the 
electron self-energy, i.e. to impose the Born Approxi- 
mation (BA). Within the first Born Approximation, the 
self-energies are evaluated with the unperturbed Green's 
functions above and obtained by the Feynman rules as 
follows in 
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huj) + D+^{uj)G^{E - huj)]M^diJii) 

We neglect here the renormalization of the phonon modes 
due to the effect of the electrons which would appear via 
a self-energy analogous to those in Eqs. (H))-®. This is 
appropriate when the mass of the ions is sufficiently large 
such that Migdal's theorem holds [l^l; however, the sub- 
sequent dispersion of the phonon Green's functions in en- 
ergy space, which leads to a finite lifetime, cannot hence 
be taken into account. To go beyond the BA one can per- 
form a self-consistent procedure for the electronic Green's 
functions such that the Green's function which satisfies 
the Dyson and Kcldysh equations and that used to evalu- 
ate Eqs. (IZ])-® are equivalent. This procedure is known 
as the self-consistent Born Approximation (SCBA). 
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The self-consistent Green's functions thus obtained 
may be used to calculate properties of interest, such 
as the steady-state inelastic current and the power in- 
jected from the electrons into the vibrational modes of 
the system[i,i,[l3,|23]. 



III. CORRELATED ELECTRON-ION 
DYNAMICS (CEID) 

A. Formulation 

One of the advantages of using dynamical methods as 
a basis for electronic transport calculations is that the 
interplay between electrical properties and atomic mo- 
tion can be addressed. Conventional Born-Oppcnhcimcr 
molecular dynamics simulations enable the calculation of 
current-induced corrections to atomic forces. However, 
in such simulations the scattering of electrons from ions 
is purely elastic and the electronic structure for a given 
ionic geometry remains in a steady state. 

A principal goal, therefore, of extending molecular dy- 
namics beyond the adiabatic approximation is to under- 
stand phenomena in which both the electrons and ions 
depart from equilibrium, the subsequent interactions be- 
tween them, and the exchange of energy between the 
two subsystems. Correlated electron-ion dynamics con- 
stitutes an attempt to introduce correlated electron-ion 
fluctuations as low-order corrections to Ehrenfest dy- 
namics, by expanding the electron-ion quantum Liouvillc 
equation in powers of such fluctuations [H, The 
method enables the description of the energy exchange 
between electrons and ions in a non-equilibrium environ- 
ment and the dynamical response of the electron gas to 
the variations in the ionic distribution. 

The idea of the method is best illustrated by applying 
it to a Hamiltonian of the form ([T|) (now expressed in 
first quantization) , in which we consider electrons linearly 
coupled to a single harmonic oscillator, 

= + ^ + 1x1^0 (i?-i?o)2 (10) 

where Hq'^"^ = Hc^°\Ro) is the A'^n-electron Hamilto- 
nian in the presence of a classical oscillator centered 
at the equilibrium position Rq and F^^") denotes the 
electron-ion coupling operator, P, R are respectively the 
ionic momentum and position operators, and K^'~' is the 
Born-Oppenheimer spring constant of the harmonic os- 
cillator. The combined electron-ion density matrix pei 
satisfies the quantum Liouvillc equation: 

Pel ~ -r^[Hci, Pol]- 



Tracing over ionic degrees of freedom leads to the fol- 
lowing set of coupled equations of motion: 
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where R = R{t) is the classical mean trajectory of the 
oscillator. pi^°^ — Tri{pci} is the A'e-particle elec- 
tronic density matrix, /i'^^"'' = Tri{Ai?jOei}, A^''^"^ = 
Tri{APpei}, AR = R^ R, AP ^ P - P, AP(^-) = 
F(JV„) _ Tr{pr"^P(^")}^, F = P ^ Tr{4'^"'F(^'^)} - 
K^°{R-Ro), P = MR. 

In order to obtain a closed set of equations, the right- 
hand sides above are truncated to lowest non-trivial order 
in the clcctron-ion coupling. Thus, in the second term in 
p2p and in the second term in p3p . we make the mean- 
field approximations 

^) = Tri{(Ai?)2p,i} « ((AP)2)pW.) ^ C««p(^ei55) 

xf°^ = iTri{{Ai?,AP}pei}«i({Ai?,AP})p(^^) 

= CRPpi''^l (16) 

The above equations are many-electron equations of mo- 
tion; these are reduced to one-electron form by tracing 
out all other electrons with the help of a Hartree-Fock 
approximation to the two-electron density matrix. This 
procedure is described in detail elsewhere [l3| and leads 
to the one-electron equations of motion 

;5e = Uhc{R),Pc]-^[F,p], (17) 
in in 

h = l[i7^(i?),^]_lCfl^[F,pJ + A, (18) 
in in M 

X = ^[H,{R),\] + liFp, + pJ)-p,Fp, 
in 2 

~^Cnp[F,p,]-K^Ofi, (19) 

where all operators are now one-electron operators, and 
where we define 

Pc= N,TT,,2...Nji''''\ 

p = N,Tt,^2...nJ^''^\ 

A==iVeTre,2...7V„A(^"^ 

B. Weak scattering limit 

In this section, we make an approximate, although 
revealing, connection between the steady-state limit of 
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the CEID equations above and the SCBA. We assume 
that the vibration is in an oscillator eigenstate with A^p^ 
phonons, where K^'^Crr = {Np\i + ^)MlQ and Crp ~ 0. 
We imagine that the phonon-free electron system has set- 
tled in a steady-state with a one-electron density matrix 
Pe = / h{E)dE, where p,{E) = E« l«>/a'5(-B - EM 
is the energy-resolved density matrix. For an infinite 
open current-carrying system, the one-electron states \a) 



with energies Ea would be Lippmann-Schwinger scatter- 
ing wavefunctions with occupancies /q, set by the battery 
terminals. We ignore variations in R relative to the equi- 
librium position i?o (hence Hc{R) = Hq above), and we 
can then solve (fTSjl - lfTO)) to lowest order in F. This is 
done in [13] ■ Taking the long-time limit of the result for 
fi gives 



-F, 



ph 



-t- 



ph 



1 



Ea — Ej3 + ftf2o — ie Ea — Ep — hQo — ie 



E„ — E, 



■0 ■ 



(20) 



where e — ^ 0+ and Qq = K^'-' /M. This expression may where Gq {E) is the advanced phonon-free electronic 

be substituted into Eq ([T7)) . The commutator [F,fi\, Green's function and the self-energies Epj^^(i?) are given 

which describes the electron-phonon coupling, then be- by 
comes: 

= j[p,{E)t-^{E)-h.c.]AE 

~hi J Pph(^^)Go (i?) - h.c.]di?, (21) 



t±(E) ^ ^ VFlgA + (^^Ph + l)/a ^ph(l-/o) iVph/a V 1^ 

P""^ ' 2Af 17o ^ ' ' \E ~ Ea- hno±ie E - Ea + Mlo±ie E~Ea + hno±ie E - Ea - hQa ± ie J ^ ' ' 

±<^{E) = 27ri^^J2F\a)[{Npi, + l)5{E + hno-Ea)+NpiAE-hno-Ea)]fa{a\F. 

(22) 



It is shown in Appendix |^ that these expressions for 
the self-energies are the same as those in the first Born 
approximation. We have thus established that in the 
limit of weak electron-ion coupling, the CEID and SCBA 
steady states agree. 



To derive this we now examine the following elastic scat- 
tering problem. We consider non-interacting electrons 
coupled linearly to an infinitely heavy classical degree of 
freedom X, with some time- independent statistical dis- 
tribution x(^). We have the one-electron Hamiltonian 



C. Large Mass Limit in CEID 



In this section, we examine the limit of infinite mass, 
in which Eqs. (jlTHlOp can again be solved analytically. 
In that case R = R{0) = i?o is a constant, the equations 
of motion for fi and A decouple, and (jlTlllSp reduce to 



H{X) = Ho- FX. 



(24) 



ihpe 

ihjj. 



Imagine solving the Liouville equation ihp(X, t) = 
(23) [H{X), p{X,t)] for the one-electron density matrix 
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piX,t). Define 

p,it) = I p{X,t)x{X)dX 



Then, 



m = J Xp{X,t)x{X)dX 
A2W = I X^p{X,t)x{X)dX. 



ihpc = [Ho,pc] - [F,fl] 
ihp = [-ffo, A] - [^>M2]- 



(25) 
(26) 
(27) 



(28) 



Now consider the distribution 



This is necessary in order to take account of the screening 
of the bare ion-ion interaction by the electron-ion interac- 
tion and the corresponding contributions to the effective 
stiffness. Since the CEID equations are derived from the 
bare interaction potentials in the system, the effective 
stiffnesses, phonon modes and frequencies are no longer 
an input, but are generated as part of the simulation. 
Furthermore, the inclusion of second-order electron-ion 
coupling (via K) arises naturally from the second-order 
expansion. A reformulation of the CEID expansion for 
systems with strong electron-nuclear correlations is de- 
veloped in [2l|- The full set of one-electron equations 
of motion, including equations of motion for the ionic 
variables R, P,Crr,Crp, and Cpp = {(AP)^), are re- 
produced in Appendix [BJ 



x{X) = -[SiX-a) + 6iX + a)]. 



(29) 



We then have exactly 



CRRPe; Crr = j X^x{X)dX = a^. 



and equations ((28)) reduce identically to ((23)) . Therefore, 
in the large-mass limit, CEID is algebraically equivalent 
to the elastic scattering problem defined by Eqs. (|24p and 
((M)) . This equivalence will be used later to benchmark 
the approximate OB method used in CEID. 



D. General CEID Equations 

The original formulation of CEID, which will be used 
for the calculations in Section IIV) starts from a more 
general Hamiltonian than that in equation ((TU)) . We start 
formally from the full electron- nuclear Hamiltonian iJoi, 
which we partition as in Section [III Al H^i = Hi'^'^ (R) + 
Ti + Hi{R), where Hc^'\r) includes the bare electron- 
ion interaction, Ti is the nuclear kinetic energy operator, 
and Hi is the bare ion-ion interaction potential. Within 
the weak-coupling approximation considered here, this 
Hamiltonian is expanded about the mean ionic trajectory 
R to second order in AR. 

H,i « H^^^^R) + Hj{R) ~ (F^^^^R) + Fi{Rj) ■ AR 



-lik^^'^HR) + Ki{R)){ARj' 



(30) 



where K^^-^^R) = d^Hi^°\R) / dR^ , Fi{R) = 
-dHi{R)/dR, and Ki{R) ^ d^Hi{R)/dR^. This Hamil- 
tonian is inserted into the full quantum Liouville equa- 
tion, and an analogous procedure to that which led to 
Eqs. (|17m9p is undertaken. The reduction of the equa- 
tions of motion to one-electron form requires an ex- 
tension to the Hartree-Fock approximation to the two- 
electron density matrix to allow for the essential 
non-idempotency introduced by electron-ion correlations. 



E. CEID with Open Boundaries (OB) 

The CEID calculations below use the open-boundary 
method described in We consider a finite, though 

possibly large, system S = LCR consisting of electrodes 
L and R with a region C between them. All dynami- 
cal scattering is assumed to be confined to C. Each fi- 
nite electrode is embedded in, and weakly coupled to, 
a sea of external probes P. Probes coupled to L{R) 
are maintained at electrochemical potential p-LiR), with 
corresponding Fermi-Dirac distributions fLi^R){E). The 
open-boundary equations of motion for the one-electron 
operators p^, p., X in S are 



ift(j= +A(') (7-pe,A,A, 



(31) 



where A^"?^ denotes the electron-ion dynamical scatter- 
ing terms, and denotes the open-boundary driving 
terms. These driving terms are 



[f:<{E)G-{E) - G+{E)±<{E)]dEi32) 



where 



= Ti^lLTi^ifl, 

S<(i?) = ^fL{E)iL + ^fR{E)h 

G^iE) = {E-Ho-±'^±iA)-^, 



(33) 
(34) 



(35) 

(36) 
(37) 



where Hq = Ho{Ro) is the phonon- free one-electron 
Hamiltonian and 1 m denotes the identity operator in re- 
gion Af . 

These equations are obtained by making two approxi- 
mations. The first is to take the wide-band limit in the 
external probes P. This makes the SP coupling strength, 
r, an energy-independent parameter and the extraction 
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terms (the first two terms in Eq. (|32p ') temporally local. 
The second approximation is the introduction of a de- 
phasing mechanism in the SP coupling, characterized by 
an energy scale, A, and a dephasing time ta = h/A. Pro- 
vided C is long enough, so that ta is less than the time 
for signals to travel between L{R) and C, the dephasing 
mechanism breaks the coherence between injection into 
L{R) and subsequent scattering in C. This in turn has 
the effect of making the injection terms (the second two 
terms in Eq. ([5^ ) independent of the dynamical scatter- 
ing in C. Otherwise, the Green's function in the injection 
terms would contain a self-energy describing the scatter- 
ing in C. 

The resultant open-boundary scheme has the benefit 
of being temporally local. However, the cost is that the 
dephasing mechanism above is in turn equivalent to re- 
placing the true Fcrmi-Dirac distributions in the probes 
P by effective distributions with an energy broadening 
^ 2A, resulting in a corresponding loss of energy resolu- 
tion. The longer the device C, the smaller the broaden- 
ing, required to mask the dynamical scattering in C. In 
the calculations presented here, we have set A = 0. The 
resulting injection terms differ from those generated by 
the value of A, appropriate for a given device length, by 
an energy uncertainty that itself disappears with A. In 
the absence of phonons, A = generates the exact un- 
broadened elastic steady-state solution for the multiple 
probe battery, which in turn gives arbitrarily close ap- 
proximations to the conventional two-terminal Landauer 
steady state [ist . 

The OB method is tested by applying it to CEID in 
the large- mass limit considered in Section [III CI We take 
the electronic system to be a resonant trimcr, described 
in more detail in the following section, within a Is tight- 
binding model with non-interacting electrons. The re- 
sults obtained from the CEID calculations are compared 
to the exact static elastic steady-state, which can be cal- 
culated separately, within numerical precision, from the 
Landauer formalism, and which, in the absence of any 
approximation in the OBs, must agree identically with 
the large-mass CEID steady-state. The two steady-state 
currents as a function of effective cross-section Crr are 
presented in Fig. [l]for a variety of biases, with excellent 
agreement. 



IV. RESULTS 

An electron within a resonant molecule characterised 
by an energy width SE will have a lifetime t ^ h/6E. If 
this lifetime is sufficiently large, the electron may be ex- 
pected to undergo several electron-phonon interactions, 
which may lead to high excitation of the vibrational 
modes of the molecule. Such multiple electron-phonon 
scattering events lie beyond lowest-order perturbation 
theory. However, since the SCBA effectively sums the 
low-order scattering events to infinite order, it will cap- 
ture at least some of the pertinent physics. In view of the 
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FIG. 1; Steady-state currents as a function of efTective cross 
section for a single degree of freedom of infinite mass, for a va- 
riety of biases, in the CEID approach. The device length used 
here was 401 atoms, with 100 atoms in each electrode, with 
the parameters F and A set to 0.4 eV and 0.0 eV respectively. 
The point at which the current drops to zero corresponds to 
one of the hopping integrals in the Hamiltonian (|24p going to 
zero. (Color online) 




FIG. 2: Linear trimer weakly coupled to two one-dimensional 
electrodes. In the simulations considered here, only the cen- 
tral atom of the trimer is allowed to undergo vibrational mo- 
tion, ti is the nearest-neighbour hopping matrix element in 
the metal electrode, t2 is the electrode-trimer hopping inte- 
gral, and is is the intra-trimer hopping integral. (Color on- 
line) 

equivalence of CEID and the SCBA in the weak-coupling 
limit established in the previous section, we conjecture 
that the CEID equations will be able to capture the phe- 
nomena of interest. Neither calculation can be expected 
to be correct in the limit of strong electron-phonon cou- 
pling, or in the limit of high phonon excitation (in which 
case, effects of anharmonicity would also be significant). 

In this section, we compare SCBA and CEID for the 
following model resonant system. We consider a linear 
trimer, illustrated in Fig. [21 which is weakly coupled to 
two one-dimensional perfect metal electrodes. The cen- 
tral atom of the trimer is treated as a dynamical quan- 
tum ion, allowed to move longitudinally. We assume non- 
interacting electrons throughout. 

Due to the absence of electron-electron screening, only 
electron-ion interactions are present in the Hamiltonian; 
these are described via a single-orbital tight-binding 
model. In the CEID simulations, the ion-ion interactions 
(which are required to calculate the dynamical matrix) 
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are described through a repulsive pair potential, with 
both sets of parameters fitted to bulk gold The 
bond length in the electrode is 2.5 A which corresponds 
to a hopping matrix elements of ti ^ —3.88 eV. The 
electrode-trimer distance is 3.509 A, corresponding to a 
hopping integral of t2 ^ —1.00 eV, while the intra-trimer 
bond length is also 2.5 A (i.e., ~ —3.88 eV). All onsite 
energies arc set to zero. The electron-phonon coupling 
matrix M u sed in th e SCBA calculations is chosen to 
be M = -^/h/2Am^F{Ro), derived from the same TB 
model as that used in CEID. The CEID calculations have 
been carried out with our parallel computer code pDI- 
NAMO nil, an implementation of the CEID formahsm 
developed to run on massively parallel computers. 

With the present parameters and a band-filling of 0.5, 
a resonance of width '--^ 0.54 eV centred at the Fermi en- 
ergy appears in the clastic transmission function. Based 
on considerations of the electron Fermi velocity and the 
geometry of the resonance, for the ionic mass considered 
below, we expect multiple clcctron-phonon interactions 
in the time interval corresponding to this width. 



A. Inelastic I — V characteristics in the externally 
damped limit 

The first comparison between the two methods is to 
calculate the low-temperature inelastic correction to the 
current-voltage spectrum for the trimer in the exter- 
nally damped limit. We assign a mass of 1 atomic mass 
unit (amu) to the moving atom, such that its Born- 
Oppcnheimer vibrational frequency Hq is hfto ^ 0.20 
eV. In the OB CEID calculations, the total number of 
atoms in the chain is 601, with 100 assigned to each elec- 
trode, with the probe-electrode coupling F = 0.4 eV. The 
second-order variables C^jj , Cpp arc set to those of the 
vibrational ground state, and kept "frozen" throughout 
the simulation. This, therefore, corresponds to the limit 
of perfect dissipation of energy away from the phonon 
modes. In the SCBA calculations, the occupation num- 
ber of the phonon mode was effectively kept at A'ph ~ 0. 

The current- voltage spectra obtained for the two meth- 
ods are shown in Fig.[3l together with the second deriva- 
tive of the inelastic contribution to the current. Both 
methods capture the inelastic feature at the correct fre- 
quency and the overall drop in the conductance is sim- 
ilar. The feature obtained using the CEID calculations 
is rather broad; this results from the absence of an ef- 
fective phonon contribution to the electron self-energy in 
the OB formalism, and from finite size effects (since the 
energy levels of the system are discrete). The width of 
the SCBA feature is a result of finite electron tempera- 
ture as well as the numerical procedure for obtaining the 
second derivative. 
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FIG. 3: Current-voltage spectrum and second derivative of in- 
elastic current for the trimer system in the externally damped 
limit. (Color online) 

B. Heating and Equilibration 

The trimer system explored here exhibits two notable 
characteristics with regard to the heating of its vibra- 
tional mode. The effective phonon occupancy obtained 
at high voltages {V ^ hilo) is significantly higher than 
that obtained in a ballistic atomic chain. In addition, 
the time taken for the vibrations to equilibrate with the 
electron gas is long [2^1. The origin of these properties 
lies in the resonant character of the system, and can be 
understood, at least qualitatively, within the FGR (see 
Appendix [C| . 

To simulate heating in the CEID calculations, we "un- 
freeze" the phonon modes and allow the vibrational de- 
grees of freedom to respond to the current-carrying elec- 
tronic structure. We take the total vibrational energy to 
be Cpp/M, which assumes equipartitioning of the vibra- 
tional energy between the kinetic and potential energies 
and includes the zero-point energy [2^. Fig. |4] shows the 
current and total vibrational energy as a function of time 
for various applied voltages. At voltages below the inelas- 
tic threshold the total vibrational energy remains flat. 
Above the inelastic threshold, as the voltage increases, 
we see the two features of the heating mentioned above; 
the total energy increases greatly and the time taken for 
equilibration also increases. Furthermore, at high volt- 
ages, the current traces "cross over" , an indication of the 
onset of negative differential resistance. 

In order to compare the results here with the SCBA, 
we extrapolate the vibrational energy as a function of 
voltage to infinite times, and compare those with the 
maximal vibrational energy obtained in the SCBA, in the 
undamped limit. These results are illustrated in Fig. [51 
together with the maximal heating according to the FGR 
for the present system, and for a quantum ion of the 
same frequency in a ballistic chain. It is seen that the 
SCBA and CEID are in good agreement up to very high 
voltages; for the highest voltage on the plot, the cffec- 
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FIG. 4: Total vibrational energy of a single dynamical ion 
as a function of time and corresponding electronic current, 
for various electrochemical potential differences in the CEID 
approach. (Color online) 
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FIG. 5: Maximal vibrational energy as a function of voltage 
for the resonant trimer system, in the dynamical CEID calcu- 
lation, the SCBA and lowest-order perturbation theory. Also 
shown is the vibrational energy obtained within the FGR for 
a single ion (with the same angular frequency) in a ballistic 
chain. (Color online) 



tive phonon occupancy is iVph ^ 10. In both cases, the 
maximal vibrational energy significantly exceeds that for 
the ballistic chain. The FGR calculation on the present 
system deviates from the other methods at high bias, in- 
dicating that CEID and the SCBA are capturing higher- 
order processes that arc absent from the lowest-order 
perturbative treatment. There are two regions of dis- 
agreement between CEID and the SCBA; at voltages just 
above the inelastic threshold the heating obtained from 
the CEID calculation is lower than that in the SCBA 
or FGR calculations. We speculate that this is due at 
least in part to the width of the inelastic spectral fea- 
ture introduced by the inexact OB method used in the 
CEID calculation; the full effect of the oscillator is grad- 
ually seen with increasing bias. Secondly, at very high 
voltages, the increase in vibrational energy in the CEID 
calculations tapers off, which may be due to the explicit 
inclusion of second-order electron-ion coupling. 

As a further comparison of CEID and the SCBA, we 
can consider the asymptotic values of the inelastic cur- 
rents in the maximally-heated limit. These arc presented 
in Fig. [6] and again there is good qualitative agreement 
between the two methods. In particular, both methods 
demonstrate that negative differential resistance will oc- 
cur in this system in the undamped limit, although the 
methods predict a slightly diffferent value for the voltage 
at which the maximum current is achieved. 



C. Inelastic current as a function of cross-section 
and ionic mass 

We now turn our attention to making a direct com- 
parison of the inelastic scattering rates produced by the 
two methods. We consider externally damped conditions 



" 20 - 




Voltage [eV] 

FIG. 6: Asymptotic values of the current as a function of volt- 
age for the resonant trimer system, in the maximally-heated 
limit, for the CEID and SCBA methods. Also illustrated are 
the steady-state currents in the fully damped limit, illustrat- 
ing the effect of the heating. (Color online) 



and assign the same fixed value of to each calcula- 
tion corresponding to an oscillator cigenstate with A^ph 
quanta. The steady-state current for a given bias (IV) is 
calculated as a function of iVph/\/M and we additionally 
examine how these currents vary over a range of masses. 
The results are shown in Fig. [T] We can see that the two 
methods remain in close agreement all the way to the 
point where the inelastic current has been suppressed 
by more than 50% relative to its value for the vibra- 
tional ground state (iVph = 0). The ionic vibrational 
energy where more significant disagreements appear (for 
A^ph/VM - 10) is of the order of 2 eV. 
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FIG. 7: Inelastic current as a function of Nph_/\/M for a bias 
of IV for a variety of ionic masses. (Color online) 



V. CONCLUSIONS 

In this paper, the first direct analytical and numer- 
ical comparison has been made between the correlated 
electron-ion dynamics and the Self-Consistent Born Ap- 
proximation. The formalisms were reviewed, and clear 
connections between the schemes were presented. From 
a numerical point of view, a model system of a linear 
trimer weakly coupled to two electrodes was studied, and 
the results are in good agreement over a range of condi- 
tions, indicating that both methods describe the under- 
lying physics in a similar manner. 

Differences begin to emerge in the limit of high thermal 
excitation, suggesting that, as methods for generating 
and summing an effective scattering series to infinite or- 
der, CEID and SCBA do ultimately differ. Furthermore, 
the effective Hamiltonians that the two methods use dif- 
fer. The SCBA Hamiltonian is a sum of an unperturbed 



electronic Hamiltonian, an unperturbed phonon Hamil- 
tonian, and a linear electron-phonon coupling. In CEID, 
the Hamiltonian consists of an unperturbed electronic 
Hamiltonian, a Hamiltonian for the bare nuclei/ions (not 
the phonons), and the electron- nuclear interaction. Un- 
der extreme conditions, this difference becomes exacer- 
bated. 

However, by applying CEID to the SCBA Hamiltonian, 
it was shown that CEID generates effective self-energies 
that, to lowest order in the electron-phonon coupling, re- 
turn the Born Approximation. A challenge for further 
work, therefore, is to seek a general diagrammatic for- 
mulation of CEID that can be compared with NEGF to 
higher orders. One possible application of a diagram- 
matic expansion of CEID are statically disordered media, 
in which the effects of multiple coherent scattering are 
significant. At a practical level, the agreement found be- 
tween the two methods opens up the exciting possibility 
of combining the first-principles electron-phonon Hamil- 
tonians, developed for use in the SCBA, with the CEID 
equations of motion, in order to generate corresponding 
dynamical electron-phonon simulations for molecular sys- 
tems. 

The observed behaviour of the resonant system studied 
here is limited by the absence of electron-electron corre- 
lation in the present model calculations. Nevertheless, 
on the understanding that electron-electron interactions 
(as well as vibrational coupling to the electrodes) may 
modify this behaviour, the results give the tentative in- 
dication that if, under high enough bias, the voltage win- 
dow engulfs an electronic resonance, with the quasi Fermi 
levels of the electrodes lying in regions of low DOS, then 
enhanced phonon relaxation times and local heating in 
the resonant structure may occur, with a resultant loss 
of mechanical stability. 
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APPENDIX A: BORN APPROXIMATION 
SELF-ENERGIES 

In this appendix we derive the phonon contribution 
to the electron self-energies within the first Born Ap- 
proximation, and show that they are equivalent to those 
obtained by substituting a current-carrying steady-state 
density matrix into the CEID equations. We consider 
here only a single ionic degree of freedom, with mass M, 
vibrational frequency f2oi and electron-phonon coupling 
matrix 



The phonon-free electronic Green's functions are 

G<{E) = 27ri^|a)/,5(^-S,)(a|. (A3) 

Hence, from Eq (l7|), 
t<^{E) ^ 27ri^M|a)[(7Vph + l)(5(S + fir!o-£;a) 

Q. 



which is the same as ([22)) . 

By combining the first and third term in ([5]), and by 
performing the contour integration, one obtains 



— / iD<{u;) + D+iij))G+iE - hw)du; 



(iVph + l)|a)(a| 



Npi,\a){a\ 



^E - Ea - hQo + ie E - E^ + hQo + ie 
The second term in (|5]) gives 



— / DU^)G<{E~nuj)diu 



1 



E - Ea- hQo + ie 



1 



E - Ea + hflo + \(., 
By combining these terms, we obtain 

A^ph + 1 - /a 



E- Ea-hno± ie 

A^ph + fa 



E - Ea + hflo ± ie 
which is the same as Eqs. 



a\M (A4) 



APPENDIX B: CEID EQUATIONS OF MOTION 

The one-electron CEID equations of motion for non- 
interacting electrons are: 



P 

B - — P - F 

iVlji 



(Bl) 



F, = Fl + Ty{p,F,} - E ^T{k,ypL,,} (B2) 

Pc = ^ [^^o, Pc] - ^ Et^'" 



(B3) 



1 



1 



C 



PcF^Pc + E ^RR iPiy"^T^{F^P'^'} - P'i>"FyPu') 

v' v" 

u' u' 

^{fl^> Ktyi,> Pc + PcK,y^' fljy,) (B5) 



RR 



Cp\ = ^ + Tr{F.A.'} - E ^--"C'RR 

v" 

C'p'p = Tr{i^^Ay' + Ai/i^y'} 



(B6) 
(B7) 

(B8) 



Above, R^,Pu are the mean position and momentum of 
the i^th ionic degree of freedom, of mass Mi, and 

= ^-s i^vij' = 



dR^dR^. 



The second-order ionic variables are 
C]^-;; = TreTri{(Ai?,Ai?,0/5ci} 
CiTp = TreTrii{(AP,Ai?,, + Ai?,, AP,)pei} 

= TreTri{(AP,AP,0/5oi}. 
ZJfl/f is defined as the inverse of Cpp^ such that 

^RR^RR —0^1,'. 

u" 

Finally, i?,,. = /v^,, + Tre{X,,'Pe}. 

These equations, along with the OB formalism de- 
scribed in Section UlI El are those used in Section HVl 
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APPENDIX C: HEATING WITHIN THE FERMI 
GOLDEN RULE 

As mentioned in the main text, the resonant system 
considered here exhibits enhanced heating under bias, 
with correspondingly large phonon relaxation times. In 



J 



this section, we examine these phenomena qualitatively 
within first-order prturbation theory. Within the Fermi 
Golden Rule, if the electron-phonon interaction is con- 
sidered as the perturbation, one can estimate the rate 
of energy transfer (the power) U injected into a single 
vibrational mode of angular frequency ilo 



U = 



27rft(iVph + 1) 
M 

2nhNpi, 



M 



^ / dEf^{E)[l - Jp{E - hno)]Tv[b^{E)Fbp{E - hno)F] 



a,p=L,R 



f dEU{E)[l - fp{E + hno)]Tv[b^{E)Fbp{E + hn„)F] 



(CI) 



a.,l3=L.R 



r 



where M is the mass of the ionic degree of freedom, 
and indices a, (3 label the Lippmann-Schwinger scattering 
wavef unctions, originating from the respective electrodes, 
with occupancies fa{E) in the Landauer picture. Da{E) 
is the partial density of states operator for the respective 
class of states. F is the electron-phonon coupling oper- 
ator discussed earlier. Defining U = Npi^hflQ, Eq. (jCip 
can be rewritten as 



where Taf3{E) = TT[Da{E)Fbp{E)F] and hl,r are the 
electrochemical potentials of the left and right battery 
terminals. For a system with reflection symmetry about 
the origin (which we can assume here), the terms in the 
integral cancel identically, and hence k and wq are given 
by 



U{t) = -KU{t)+W(). 



(C2) 



For the FGR calculations of maximal heating in Fig O 
the quantities k and wq were computed by full energy 
integration, from Eq. (|Cip . with the maximal heating 
being given by the zero power condition J/max = wq/k. 

For the purposes of gaining physical insight into the be- 
haviour of the resonant calculation, let us now simplify 
the calculation as follows: we assume zero electronic tem- 
perature, and assume that the variations in the electronic 
Green's functions over energies of the order of hTlo are 
small such that 

TT{bL{E)FbR{E ± hno)F} 

« Tr{£)i,(£; ± hQ.(i/2)FbR{E ± hVtn/2)F} 

±!!:^Tr{Z)L(-B ± hno/2)Fb'i^{E ± hno/2)F} 
zfl^Tiib'^iE ± hna/2)FbR{E ± hno/2)F} + o{hna 
Hence, for /ii — fiR > hilg, 

"^^^ [TLhit^h) + Timinii) + TLRiiJ-L) + Tlr{ij,r)] 



M 



[Tr{bL{E)Fb'^{E)F} 



^ [TLLifJ-L) + Timi^ii) + TLRiflL) + TLRiflR) 



2TTh /•A'i-«2o/2 



Wo 



M 



Tlr{E) dE 



(C3) 



-TT{b'L{E)FbR{E)F}]AE + Oih^o 



Consider now a situation in which under a large enough 
bias, the energy window for conduction engulfs an elec- 
tronic resonance, with /i^ and now lying in regions 
of low densities of states, on each side of the resonance. 
Then k, which collects contributions from energies in the 
vicinity of the two Fermi levels, is small, and gets smaller 
|2with increasing bias, as and ^.r move further away 
from the resonance. Since k depends quadratically on the 
density of states, its decrease with bias should be faster 
than linear, wq, on the other hand, collects contributions 
from the entire conduction window, and saturates with 
increasing bias. This is the origin of the enhancement of 
C^max and of the phonon equilibration time t = com- 
pared with the ballistic case (in which wq increases lin- 
early with bias, and k is approximately bias-independent 
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